
function [r, rt, w] = tremble_r(param, util_type, N, smooth_wind, T, yi, reset_shocks, rLower, rUpper, conv_r, counter_r, maxiter_r, crit_r, precise, stateGrid_size,eps_eta_mat,sigma_tremble, p_unemp, w_unemp)

bta = param(1);
r = param(2);
abar = param(3);
W_cc = param(4);
kappa = param(5);
deltax = param(6);
sigma_c = param(7);
psi = param(8);
cov_fun = param(9);
gam = param(10);
sigma_l = param(11);
lbar = param(12);
delta_eta = param(13);
delta_inh= param(14);
delta    = param(15);
alpha    = param(16);
b        = param(17);
btat     = param(18);
amin     = param(19);
amax     = param(20);
na       = param(21);

% keep track of interest rate and market clearing in each iteration

r_all = NaN(1,maxiter_r);
r_all(1) = r;
K_all = NaN(1,maxiter_r);
A_all = NaN(1,maxiter_r);

tic
while conv_r > crit_r && counter_r < maxiter_r
    
    %fprintf('r = %d\n', r)
    
    [conv_crit, A, K] = tremble_r_fp(r,param, util_type, N, smooth_wind, T, yi(1:T,1:N),reset_shocks(1:T,1:N), stateGrid_size,eps_eta_mat(1:T,1:N),sigma_tremble, p_unemp, w_unemp);
    
    % Step 5: compute excess demand/supply of asset and use bisection
    % method to adjust r accordingly
    if conv_crit > 0
        rnew = 0.5.*(r+rLower);
        rUpper = r;
    else
        rnew = 0.5.*(r+rUpper);
        rLower = r;
    end
    conv_r = abs(r-rnew);
    
    r = rnew;
    counter_r = counter_r + 1;
    
    r_all(counter_r+1)=r;
    A_all(counter_r)=A;
    K_all(counter_r)=K;
    
end
toc

[~ ,min_index] = min(abs(A_all-K_all)); % get the iteration with best mkt clearing
r = r_all(min_index);
A = A_all(min_index);
K = K_all(min_index);
w = (1-alpha).*K.^alpha; % Wage (from firm l-FOC)

if precise == 1
    
    
    obj = @(x) tremble_r_fp(x,param, util_type, N, smooth_wind, T, yi(1:T,:),reset_shocks(1:T,:), stateGrid_size,eps_eta_mat,sigma_tremble, p_unemp, w_unemp);

    r0 =r;
    
    
    %%%
    fsolve_options = optimoptions( 'fsolve');
    fsolve_options.Display = 'iter';
    fsolve_options.MaxIterations = 100;
    fsolve_options.StepTolerance = 1e-4;
    %fsolve_options.TypicalX = 0.1*ones(length(x0),1);
    fsolve_options.FunctionTolerance = 1e-4;
    
    
    tic
    [xout, ~] =  fsolve( obj, r0, fsolve_options);
    toc
    
    [~,A,K]  = obj(xout);
    w = (1-alpha).*K.^alpha; % Wage (from firm l-FOC)
    r = xout; 
end

rt = (r*(1-deltax*delta_inh) + deltax*(1-delta_inh))./(1-deltax); % effective interest rate (includes insurance market)

disp(['                                r        rt        A        K        w ']) 
disp(['   Trembles        :        ' , sprintf('%1.4f   ', [r, rt,A, K, w])  ]);
